Self-Organized States in Cellular Automata: Exact Solution 
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The spatial structure, fluctuations as well as all state prob- 
abilities of self-organized (steady) states of cellular automata 
can be found (almost) exactly and explicitly from their Marko- 
vian dynamics. The method is shown on an example of a 
natural sand pile model with a gradient threshold. 

PACS numbers: 05.40. +j, 03.20.+i, 46.10.+Z 



Complicated dynamics of various discrete systems may 
naturally be modeled by CA having rather simple itera- 
tion rules. In particular, CA models are useful to study 
traffic jams H, granular material dynamics and self- 
organization|3|4|] . The transport properties of systems 
the, especially in the Self-Organized (SO) critical regime 
were extensively studied (cf., However, for many 

applications, one needs to know an entire structure of 
the SO states. These may include: (i) the temperature 
profiles for the convection dominated thcrmoconduction 
and turbulent convection || (as in the convective zone of 
the sun and stars), (ii) stability and average profiles of 
granular materials (e.g., the sand pile profiles) and gran- 
ular flows [^|J^1, (iii) equilibrium and steady state pro- 
files of plasma pressure and temperature in fusion devices 
flllOj, etc. Despite its importance, the problem of spa- 
tial structure and characteristics of the SO steady states 
has received minor attention |Q,0-O|. In this Letter, we 
propose a method which is equally applicable to any CA 
provided its rules are rather simple. We illustrate it on 
the simplest and popular example of an one-dimensional 
sand pile with a gradient stability criterion. Note that 
this model is the closest to a natural pile of sand created 
by random sprinkling of sand and with a known repose 
angle. We focus ourselves on calculating of an average 
slope (the calculation is more transparent) though calcu- 
lating other characteristics is trivial. 

An interesting result obtained is that the SO profiles 
of a local slope are non-trivial even for this simplest case. 
They are typically flat (e.g., linear) throughout the pile 
while a narrow region (similar to the boundary layer) 
with rapidly increasing gradient always occur near the 
top of the pile. This picture is very similar to that ex- 
perimentally measured in a strongly turbulent convection 
of a passive scalar (temperature) with a mean gradient 
Pj . The region of the super-critical (unstable) gradient 
may form near the bottom of the pile when the noise is 
strong enough to maintain almost continuous sand flow 
(overlapped avalanches). This is in excellent agreement 
with direct numerical simulations JlO| . 

Among many sand pile models, the Abelian model is 
the only one that was proved to be analytically exactly 



solvable JlJJlJ]. It is, however, rather unnatural since 
its stability criterion is given in terms of a site height, 
not slope of the pile (which depends on neighboring sites, 
too). Another tacit assumption is one of weak noise: no 
sand is added to the pile unless at least one unstable site 
presents (i.e., until the very last avalanche is gone). Note 
that "waves of toppling" [ p^[ , which are the main point 
of the whole analysis, are well defined in the weak noise 
limit, only. We show that both these restrictions may 
be relaxed: (i) a site stability criterion may depend on 
site's neighbors and (ii) sand may be added at every time 
step thus affecting an avalanche dynamics. We also show 
that in the weak noise limit, all state probabilities can be 
calculated exactly in a closed form, while for the strong 
noise, they can be found almost exactly, i.e., with any a 
priori given accuracy. 

For a general iV-dimensional sand pile automaton, the 
procedure is following. (1) Reformulate the model in the 
representation where the stability is local and defined by 
the state of a site, alone. This can be done at expense 
of introducing a nonlocality to the toppling and noise 
rules (i.e., they may depend on states of adjacent sites 
too). (2) Consider the dynamics of a single site. Since 
toppling rules have no intrinsic memory, however, it is 
Markovian. Construct an iV-dimensional Markov hyper- 
lattice (an analog of a Markov chain in one dimension) 
with the transition probabilities defined by the CA rules. 
All the transition probabilities which depend on states 
of other sites are, for now, free parameters. Introducing 
a generating function, one can then solve the problem 
for a single site. (3) Noting that all sites are identical, 
we relate the Markov transition probabilities for different 
sites. Boundary conditions then uniquely define their 
values and, thus, the SO state of the pile. Note that a 
mean- field type closure is needed at step (3), only. 

Model — The model we consider consists of L + 1 spa- 
tial sites, numbered from x = to x = L. To each site x is 
assigned a variable h(x), the height of the site. CA rules 
are applied to the pile at each time-step. Sand grains are 
added randomly to sites with probabilities p sa nd{x) in- 
creasing the height by one. When the site x is unstable, 
i.e., if the local slope [difference of heights of two neigh- 
boring sites h(x) — h(x + 1)] exceeds some critical value 
Ah cr it{x), Nf sand grains topple onto the neighboring 
site x + 1 (local, limited model [Q). Note that Nf > 1 
corresponds to the physical situation where friction be- 
tween sand grains at rest is greater than friction of those 
in motion. Sand is expelled from the pile through the 
right end x = L. 

For the stability condition to be local, we represent 
the sand pile in "gradient space" . We assign to any site 



1 



the local difference of heights of the nearest neighbors 
Ah(x) = h(x) — h{x + 1). Then the noise and toppling 
rules become, respectively: 

noise I *h(x-l)^Ah(x-l)-l, , , 

nolse \ Ah{x) -> Ah(x) + 1; [1& > 

Ah(x - 1) -» Aft(a: - 1) + iV>, 
toppling <( A/i(x) -> A/i(x) - 2JV>, (lb) 
A/i(x + !)->■ Ah(x + 1) + iV>. 



The left end (x = 0) is now an open boundary (top of the 
pile): Ah(0) -> A/i(0) - 2iV/, A/i(l) -> Aft(l)+JV/; the 
right end (a; = L) is a closed boundary (bottom of the 
pile): Ah(L - 1) -> A/i(i - 1) + iV>, A/i(L) -> A/i(£) - 
iVy. Note that "particles" in the gradient space are not 
sand grains! They may enter or leave the system through 
the left boundary, only. Noise creates a "particle" at 
x = with the probability p sa nd{Q), as follows from Eq. 
( |la| ) for x — 0. A toppling at x = results in a sudden 
loss of Nf "particles" . 

Zero-dimensional pile — Now, we may consider one 
site x, alone. It is described by a collection of states 
representing all possible values of local gradient. Negative 
states are not allowed. These states are labeled by an 
integer variable k = Ahix), and the critical slope is Z c = 
Ah cr u(x). The states k < Z c are stable, those k > Z c 
are unstable, and the state k = Z c is marginally stable. 
We introduce the probabilities pt for a site to occupy a 
state k, i.e., to have the slope Ah = k. Due to noise and 
overturning events, the state of a site will evolve in time. 
The rules given by Eqs. (1) are independent of previous 
history of a system. Therefore, they define the evolution 
of the slope of a site to be a Markov process. The states 
arranged in increasing order of k form a Markov chain. 
Adding and toppling rules specify transition probabilities 
from one state to another on this chain. 

1. Adding sand [Eq. (la)] results in jumps by 1 right or 
left (i.e., a increase or decrease gradient). The transition 
probabilities of the process are a and /?, respectively, and 
equal: a = p san d(x + 1) [l - p sa nd(x)] , = p san d{x) [l - 
Psand(x + 1)] . Note here that adding of a sand grain in 
real space results in increase or decrease of a state (i.e., 
local slope) in gradient space. 

2. Toppling of a site [Eq. (lb)] results in a jump by 



2Nf states left (i.e., a decrease in gradient). The proba- 
bility of that process is 1. , i.e., an unstable state topple 
on the next time-step with the probability unity. 

We introduce two "nonlocal" transition probabilities. 
(1) Toppling of one of two neighboring sites results in a 
jump by Nf states right (i.e., an increase in gradient). 
The probability of this process is written as e* . (2) Top- 
pling of both two neighbors results in a jump by 2Nf 
states right. The transition probability is written as 5*. 
Both e* and 5* are simply constants here which are to be 
specified in the one-dimensional model via a mean-field- 
type closure. Generally speaking, a, /?, e*, and S* will 
depend on x and k. We later consider the case where all 



are independent of x and k (homogeneous pile with no 
local slope dependence). 

In Fig. pj, the Markov chain with all possible transi- 
tions for stable (•) and unstable (o) states with corre- 
sponding transition probabilities is shown. Noise results 
in one-step random walk of a particle on this chain. Top- 
pling of sites results in jumps by Nf and 2Nf states. 
Since the noise process is statistically independent of top- 
pling, these processes may combine with each other re- 
sulting in jumps by Nf - 1, Nf + 1, 2N f - 1, and 2N f + 1 
jumps, with the probabilities respectively proportional to 
e*(3, e*a, 8*/3, and 5*a. All other transition coefficients 
are similarly defined. 

We thus have reduced the problem of a sand pile to 
the problem of a random walk of a particle on a chain of 
states where the transition probabilities are exactly de- 
fined. For a general type of a Markov process the general 
kinetic (or master) equation is 



Pn(t) = {lnkPk(t) ~ lknPn{t)} 



k=0 



(2) 



where jkn are the transition probability coefficients from 
state n to state k. Note that the term ^ n kPk describes 
transitions into the state n from state fc, while "fknPn 
corresponds to transition out of n into other states k. 
This equation defines the probabilities p n for the system 
to be in state n. The general kinetic equation for one 
site can easily be written using Fig. El Because of space 
limitations, we do not write it explicitly. 

We introduce a generating function for the probability 
distribution pk '■ 



F(c,t) = j2c k Pk (t), 



(3) 



k=0 



where £ can take values |C| < 1 for a series to converge. 
The probability distribution can be recovered from the 
generating function as 



Pk (t) = (l/kl)d k F((,t)/dC k 



IC=o ' 



(4) 



Some properties 
F(M) = 1, 



of 



C=i 



the generating function are 
(n(t)),..., where 'prime' means 

derivative. Here the first is the normalization condition: 
J^Pj — 1; and second relates F((,t) to the first moment 
(e.g., expectation value) of the probability distribution. 
Higher moments (i.e., standard deviation, etc.) are ob- 
tained from higher derivatives of F(£,t). Multiplying 
each of Eqs. @ for p\ respectively by ( k and taking sum 
over all < k < oo, we straightforwardly obtain an equa- 
tion for a generating function F = + Since we are 
interesting only in a steady state, we set F = 0. Then it 
reads 

nOW + *(C) (C 2Nf i) t^ 1 + {•}] = PoiC 1 i) , 

(5) 
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where a = (3, e = e*/a, 5 = 6* /a and {•} = (( + ( 1 - 
2)(1 + ea( N f + Sa( 2N f ) + e(( N f - 1) + 5{( 2N f - 1). Here 

*(C) = E&oC*Pfe, *(C) = Etz.+iCVfe are the par- 
tial generating functions of stable and unstable states, 
respectively. 

To find one may use a simple trick [using Eq. da)]: 



*(C) 



0. 



(6a) 



In general, this system is an infinite set of related equa- 
tions. It relates all pz c +i+j to p a . Together with 
F(l) = 1, it provides an exact solution, F((), of Eq. (||). 
For an Abelian sand pile, the transition probabilities of 
simultaneous toppling and noise [e.g., toppling to higher 
states, 1.5*a (see Fig. [j])] vanish identically by defini- 
tion. (Note, this is the limit when noise is too weak to 
affect avalanche dynamics.) Thus, the highest achievable 
state is Z c + 2Nf. The Markov chain is finite, and Eqs. 
( |6a| ) reduce (schematically) to the system with triangular 
matrix (with ay and being constants) 



/ fti.i ai, 2 ■ ■ ■ ai,2N f \ 
a 2j 2 ■ ■ ■ a 2 .2N f 
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(6b) 



which can be solved exactly and explicitly. In the opposite 
case, when noise is not weak, one may, however, truncate 
the system of Eqs. ( |6a|) to a finite hierarchy simply by 
noticing that the probabilities pz a +2N f +j, j > 1 are very 
low since they can be reached only from unstable states. 
If one truncates at the state Z c + 2Nf + j, the error in 
determination of all the state probabilities, pk, will not 
exceed the value (1.5* a)? . 

Eqs. (|6a| ) or (g^) constitute the (almost) exact solu- 
tion, i.e., state probabilities of a site x, of a CA sand pile 
model in terms of state probabilities (entering through e 
and 5) of its neighbours. 

The normalization condition, F(l) = 1, gives the rela- 
tion for the total probability for a site to be unstable: 



a[p +N f (e + 2S)]/2N f . 



(7) 



We straightforwardly define the SO profile as the math- 
ematical expectation value (average) of the random pro- 
cess: 



(n}=FUl) 



(8) 



where two unknowns p${e,8) and ^(l) appear and are 
to be found from Eqs. ([3a]) or (6b). For arbitrary Z c and 
Nf, (especially for large Z c and Nf, i.e., in the contin- 
uous limit) the result can be easily found numerically. 
To obtain an analytically tractable expression, we make 
additional (though natural) approximations. 



1. Let's consider an asymmetric random walk on a fi- 
nite chain with transition probabilities to the right and to 
the left, g and r, respectively. One can easily show [recur- 

sively from Eq. §)} that p (f) = (f - l) [(f) C - 
where c is a constant which is found from an expan- 
sion of Po{g/r) near g/r ~ 1 to give Po| ff/ v=i = V c - 
By analogy with an asymmetric random walk, we write 
g/r = {a + N f e* + 2N f 5*)//3 =1 + N f e + 2N f S and 



Po 



N f (e + 26)/([l + N f (e + 2S)} 1 /Po° 



(9) 



where p^ 0) = p \e=s=o = (Z c - Nf + 3/2)" 1 [the last 
follows from Eq. (||) for e = d~ = 0], 

2. To define ^(l), we consider two limits for which 
^(C) is known a propri. When e = 5 = 0, only one- 
step transitions (noise) exist. Therefore, from the def- 
inition, we have ^(l)|e=«=o = [Z c + l)p Zc +i, i.e., only 
the first unstable state can be achieved. For e, 5 suf- 
ficiently large, the states k € [Z c + 1, Z c + 2Nf] are 
roughly uniformly populated while higher states Z c + 
2Nf + k,k > 1 have low probability, as they can be 
reached only from unstable states. Thus, we can write 
1^(1) ~ (Z c + 1 + (2N f - l)/2). From comparison 
of the last two equations, we conclude 



9' c (i) ~ V [(Z e + 1) + (N f - 1/2) f(e,5)] 



(10) 



where / = 1 for large e, S and vanishes for e = S = 0. 
Because the quantity e + 25 is a "measure of asym- 
metry" of a random walk, we choose /(e, 5) = 2(e + 
25)/{l + {e + 25)). 

Finally, the SO local slope (for reads 



Ah so 



[po + N f {e 
+ (Z, 



1 



-25)} 
3JV//2 



1 



N f - 1/2 
f JV>(e + 2<S) + 
\-5N f /(e + 25). 



(11) 



Here po is given by Eq. (||). Eq. ( |Tl| ) depends on the 
noise strength a as well as on the toppling probabilities 
of adjacent sites of the pile. 

One- dimensional pile — Eq. (11) defines the average 
SO slope for every site x. The quantities e and 5 are de- 
fined by toppling probabilities of neighboring sites. Each 
site x topples with probability V — V{x). This probabil- 
ity varies from site to site. In the mean-field approxima- 
tion, by definition 



ea =V{x-l)[l~P(x+l)] 

+V(x + l)[l-V{x-l)], 
5a = V(x-l)V{x + l). 



(12) 



Note, the stonger the noise, the better this anzatz works, 
because of decorrelation of V(x ± 1) caused by noise. Eq. 
(0) can be written as a recurrence equation for probabil- 
ities V{x) 

2V(x) = ap (x)/Nf + [P(x - 1) + V(x + 1)] , (13) 
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where po{x) is also a function V and given by Eq. (||). 
This equation can be solved numerically with the condi- 
tion at the open (left) boundary that "influx" = "outflux" . 
The initial value is thus, V(0) = p 3an d(0)/Nf. Eq. |y]) 
together with Eqs. jf^ ) defines a spatial profile of the SO 
slope of the sand pile. In the continuous limit (vanish- 
ing cell size of a Markov chain), Eq. |l^) is equivalent to 
V" x = apo/Nf. Thus the approximate solution matching 
the boundary condition is 

V{X) ~ (psandPo/Nf) X 2 + ( Psan d(l ~ Po)/N f ) X. (14) 

The SO gradient profiles are shown in Fig. || for Z c = 
8,Nf = 3, and three values of noise strength a ~ p sa nd'- 
a = 1/5000 (low noise), a = 1/1500 and a = 1/500 
(high noise). The average gradient profiles always have 
a region of relatively small gradient, "boundary layer", 
near the top of the pile. This region appears due to 
the effect of the open (in gradient space) boundary at 
x = 0. In the case of low noise, the SO gradient is al- 
most constant throughout the pile and always below the 
marginally stable value. In the "over-driven" regime (i.e., 
high noise), a region of a super- critical gradient appears 
near the bottom of the pile. In this case the sand influx 
is so large that a nearly constant flow of sand forms near 
the bottom, thus maintaining an unstable, super-critical 
gradient. These results agree well with simulations [ jTo| . 

In this Letter, we show that Abelian property is not 
necessary for an (almost) exact solvability of a sand 
pile CA. As an example, a spatial profile of an one- 
dimensional sand pile is calculated. 

We are grateful to M.N. Rosenbluth for numerous in- 
teresting and fruitful discussions and critical comments. 
We also thanks B.A. Carreras, D.E. Newman, and T.S. 
Hahm for discussions. This work was supported by De- 
partment of Energy grant No. DE-FG03-88-ER53275. 



(1989); Phys. Rev. A 45, 7002 (1992); D. Dhar and R. 

Ramaswamy, Phys. Rev. Lett. 63, 1659 (1989); H.M. 

Jaeger, et. al, Phys. Rev. Lett. 62, 40 (1989). 
[7] A. Vespignani, et. al, Phys. Rev. E 51, 1711 (1995). 
[8] PH. Diamond,T.S. Hahm, Phys. Plasmas 2, 3640 (1995). 
[9] J. P. Gollub, et. al, Phys. Rev. Lett. 67, 3507 (1991); BR. 

Lane, et. al, Phys. Fluids A 5, 2255 (1993); B. Castaing, 

et. al, J. Fluid Mech. 204, 1 (1989). 
[10] D.E. Newman, et. al, Phys. Lett. A 218, 58 (1996); Phys. 

Plasmas, 3, 1858 (1996). 
[11] J.M. Carlson, et. al, Phys. Rev. Lett. 65, 2547 (1990). 
[12] D. Dhar, Phys. Rev. Lett. 64, 1613 (1990); S.M. Majum- 

dar and D. Dhar, Physica A 185, 129 (1992). 
[13] V.B. Priezzhev, J. Stat. Phys. 74, 955 (1994); E.V. 

Ivashkevich, J. Phys. A 27, 3643 (1994). 
[14] E.V. Ivashkevich, et. al, J. Phys. A 27, L585 (1994); E.V. 

Ivashkevich, Phys. Rev. Lett. 76, 3368 (1996). 



i 

stable i unstable 




FIG. 1. The Markov chain representing a collection of 
states (slopes) for any site x in the sand pile model. 
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FIG. 2. The SO profiles of gradient of the pile 
(Z c = 8, Nf = 3) for three noise levels. 
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